#include "delaunay.h"
#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <iostream>
#include <unordered_map>
#include <deque>
#include <vector>

Point::Point() : x(0), y(0) {}
Point::Point(double x, double y) : x(x), y(y) {}

const double EPSILON = 1.0 / 1048576.0;  // 与JS版本保持一致的精度

struct Circumcircle {
    int i, j, k;
    double x, y, r;
    Circumcircle() : i(0), j(0), k(0), x(0), y(0), r(0) {}
    Circumcircle(int i, int j, int k, double x, double y, double r) : i(i), j(j), k(k), x(x), y(y), r(r) {};
};

// 计算包含所有点的超级三角形
std::vector<Point> supertriangle(const std::vector<Point>& vertices) {
    double xmin = std::numeric_limits<double>::infinity();
    double ymin = std::numeric_limits<double>::infinity();
    double xmax = -std::numeric_limits<double>::infinity();
    double ymax = -std::numeric_limits<double>::infinity();
    
    for (const auto& vertex : vertices) {
        if (vertex.x < xmin) xmin = vertex.x;
        if (vertex.x > xmax) xmax = vertex.x;
        if (vertex.y < ymin) ymin = vertex.y;
        if (vertex.y > ymax) ymax = vertex.y;
    }
    
    double dx = xmax - xmin;
    double dy = ymax - ymin;
    double dmax = std::max(dx, dy);
    double xmid = xmin + dx * 0.5;
    double ymid = ymin + dy * 0.5;
    
    return {
        Point(xmid - 20 * dmax, ymid - dmax),
        Point(xmid, ymid + 20 * dmax),
        Point(xmid + 20 * dmax, ymid - dmax)
    };
}

// 计算三角形的外接圆
Circumcircle circumcircle(const std::vector<Point>& vertices, int i, int j, int k) {
    double x1 = vertices[i].x;
    double y1 = vertices[i].y;
    double x2 = vertices[j].x;
    double y2 = vertices[j].y;
    double x3 = vertices[k].x;
    double y3 = vertices[k].y;
    
    double fabsy1y2 = std::abs(y1 - y2);
    double fabsy2y3 = std::abs(y2 - y3);
    
    double xc, yc, m1, m2, mx1, mx2, my1, my2, dx, dy;
    
    /* 检查重合点 */
    if (fabsy1y2 < EPSILON && fabsy2y3 < EPSILON)
        throw std::runtime_error("Eek! Coincident points!");
    
    if (fabsy1y2 < EPSILON) {
        m2 = -((x3 - x2) / (y3 - y2));
        mx2 = (x2 + x3) / 2.0;
        my2 = (y2 + y3) / 2.0;
        xc = (x2 + x1) / 2.0;
        yc = m2 * (xc - mx2) + my2;
    }
    
    else if (fabsy2y3 < EPSILON) {
        m1 = -((x2 - x1) / (y2 - y1));
        mx1 = (x1 + x2) / 2.0;
        my1 = (y1 + y2) / 2.0;
        xc = (x3 + x2) / 2.0;
        yc = m1 * (xc - mx1) + my1;
    }
    
    else {
        m1 = -((x2 - x1) / (y2 - y1));
        m2 = -((x3 - x2) / (y3 - y2));
        mx1 = (x1 + x2) / 2.0;
        mx2 = (x2 + x3) / 2.0;
        my1 = (y1 + y2) / 2.0;
        my2 = (y2 + y3) / 2.0;
        xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
        yc = (fabsy1y2 > fabsy2y3) ?
            m1 * (xc - mx1) + my1 :
            m2 * (xc - mx2) + my2;
    }
    
    dx = x2 - xc;
    dy = y2 - yc;
    return Circumcircle(i, j, k, xc, yc, dx * dx + dy * dy);
}

// 快速边去重优化版本
void dedup(std::vector<int>& edges) {
    const size_t edgeCount = edges.size() / 2;
    if (edgeCount == 0) return;

    // 预分配哈希表空间
    std::unordered_map<int64_t, uint8_t> edgeMap;
    edgeMap.reserve(edgeCount);

    // 计算每条边出现的次数（只需知道是1次还是多次）
    for (size_t i = 0; i < edgeCount; ++i) {
        const int a = edges[i * 2];
        const int b = edges[i * 2 + 1];
        const int64_t key = (static_cast<int64_t>(std::min(a, b)) << 32) | static_cast<uint32_t>(std::max(a, b));
        edgeMap[key]++;
    }

    // 直接在原数组上修改，避免创建新数组
    size_t writePos = 0;
    for (size_t i = 0; i < edgeCount; ++i) {
        const int a = edges[i * 2];
        const int b = edges[i * 2 + 1];
        const int64_t key = (static_cast<int64_t>(std::min(a, b)) << 32) | static_cast<uint32_t>(std::max(a, b));
        if (edgeMap[key] == 1) {
            edges[writePos++] = a;
            edges[writePos++] = b;
        }
    }

    // 调整数组大小
    edges.resize(writePos);
}

// 内存池实现，减少频繁内存分配
class MemoryPool {
private:
    std::vector<char> pool;
    size_t nextAvailable;

public:
    MemoryPool(size_t initialSize = 1024 * 1024) : pool(initialSize), nextAvailable(0) {}

    void* allocate(size_t size) {
        if (nextAvailable + size > pool.size()) {
            // 扩展内存池，按1.5倍增长
            const size_t newSize = std::max(pool.size() * 3 / 2, nextAvailable + size);
            pool.resize(newSize);
        }
        void* ptr = &pool[nextAvailable];
        nextAvailable += size;
        return ptr;
    }

    void reset() { nextAvailable = 0; }
};

// 德劳内三角剖分主函数 - 优化版本
std::vector<int> triangulate(const float* coords, size_t num_points, bool (*progress_callback)(float)) {
    std::vector<Point> vertices;
    vertices.reserve(num_points);
    for (size_t i = 0; i < num_points; ++i) {
        vertices.emplace_back(coords[2*i], coords[2*i+1]);
    }
    return triangulate(vertices);
}

std::vector<int> triangulate(
    std::vector<Point> vertices,
    std::optional<std::function<Point(const Point&)>> vertexProcessor
) {
    int n = vertices.size();
    
    // 如果点数少于3，无法构成三角形
    if (n < 3)
        return {};
    
    // 如果提供了顶点处理器，应用它
    if (vertexProcessor) {
        std::vector<Point> processedVertices;
        processedVertices.reserve(n);
        for (const auto& vertex : vertices) {
            processedVertices.push_back((*vertexProcessor)(vertex));
        }
        vertices = std::move(processedVertices);
    }
    
    // 创建索引数组并按x坐标升序排序
    std::vector<int> indices(n);
    for (int i = 0; i < n; i++)
        indices[i] = i;
    
    // 使用lambda进行排序，现代编译器可以优化
    std::sort(indices.begin(), indices.end(), [&vertices](int i, int j) {
        double diff = vertices[i].x - vertices[j].x;
        return diff != 0 ? diff < 0 : i < j;
    });
    
    // 计算超级三角形并添加到顶点列表
    auto st = supertriangle(vertices);
    vertices.push_back(st[0]);
    vertices.push_back(st[1]);
    vertices.push_back(st[2]);
    
    // 使用deque代替vector，提高插入/删除效率
    std::deque<Circumcircle> open;
    open.push_back(circumcircle(vertices, n + 0, n + 1, n + 2));
    std::vector<Circumcircle> closed;
    closed.reserve(n * 2); // 预分配空间
    
    std::vector<int> edges;
    edges.reserve(n * 6); // 预分配空间
    
    // 内存池，减少小对象分配
    MemoryPool pool;
    
    // 逐点插入（改为正序遍历，提高缓存友好性）
    for (int i = 0; i < indices.size(); ++i) {
        int c = indices[i];
        edges.clear();
        
        // 检查每个开放三角形
        for (auto it = open.begin(); it != open.end();) {
            // 如果点在三角形外接圆的右侧，则该三角形不再需要检查
            double dx = vertices[c].x - it->x;
            if (dx > 0.0 && dx * dx > it->r) {
                closed.push_back(*it);
                it = open.erase(it);
                continue;
            }
            
            // 如果点在外接圆外，则跳过
            double dy = vertices[c].y - it->y;
            if (dx * dx + dy * dy - it->r > EPSILON) {
                ++it;
                continue;
            }
            
            // 移除三角形并添加其边到边列表
            edges.push_back(it->i);
            edges.push_back(it->j);
            edges.push_back(it->j);
            edges.push_back(it->k);
            edges.push_back(it->k);
            edges.push_back(it->i);
            it = open.erase(it);
        }
        
        // 去除重复边
        dedup(edges);
        
        // 为每条边创建新三角形
        for (int j = 0; j < edges.size(); j += 2) {
            int a = edges[j];
            int b = edges[j+1];
            // 复用内存池分配的内存
            Circumcircle* cc = new (pool.allocate(sizeof(Circumcircle))) Circumcircle(circumcircle(vertices, a, b, c));
            open.push_back(*cc);
        }
    }
    
    // 将剩余的开放三角形添加到闭合列表
    closed.insert(closed.end(), open.begin(), open.end());
    
    // 构建最终的三角形列表，排除包含超级三角形顶点的三角形
    std::vector<int> result;
    result.reserve(closed.size() * 3); // 预分配空间
    for (const auto& triangle : closed) {
        if (triangle.i < n && triangle.j < n && triangle.k < n) {
            result.push_back(triangle.i);
            result.push_back(triangle.j);
            result.push_back(triangle.k);
        }
    }
    
    // 移除超级三角形的顶点
    vertices.pop_back();
    vertices.pop_back();
    vertices.pop_back();

    return result;
}

// 检查点是否在三角形内
std::optional<BarycentricCoords> contains(const std::vector<Point>& tri, const Point& p) {
    // 边界框测试，快速拒绝
    if ((p.x < tri[0].x && p.x < tri[1].x && p.x < tri[2].x) ||
        (p.x > tri[0].x && p.x > tri[1].x && p.x > tri[2].x) ||
        (p.y < tri[0].y && p.y < tri[1].y && p.y < tri[2].y) ||
        (p.y > tri[0].y && p.y > tri[1].y && p.y > tri[2].y))
        return std::nullopt;
    
    double a = tri[1].x - tri[0].x;
    double b = tri[2].x - tri[0].x;
    double c = tri[1].y - tri[0].y;
    double d = tri[2].y - tri[0].y;
    double i = a * d - b * c;
    
    // 退化三角形
    if (i == 0.0)
        return std::nullopt;
    
    double u = (d * (p.x - tri[0].x) - b * (p.y - tri[0].y)) / i;
    double v = (a * (p.y - tri[0].y) - c * (p.x - tri[0].x)) / i;
    
    // 如果在三角形外，失败
    if (u < 0.0 || v < 0.0 || (u + v) > 1.0)
        return std::nullopt;
    
    return BarycentricCoords(u, v);
}